This simulation and imaging demo illustrates the numerical and functional performance of the ASP deconvolution algorithm as implemented within CASA (currently as an experimental algorithm).
The Asp algorithm has been implemented as a minor-cycle option within the tclean task. It defines the sky model as a collection of Gaussians whose parameters are optimized as part of an iterative process. The current CASA implementation is an adaptation of the original Asp deconvolution algorithm in which each iteration optimizes/adjusts the parameters of one Gaussian and includes it in the global list of best-fit Gaussians. The algorithm begins with an automatically-generated initial guess of four scales of sizes ranging from two to eight times the PSF width, and this list is then allowed to evolve and grow. Two control parameters are offered. One is to limit the largest scale allowed for the initial guess, meant to stabilize the reconstruction for under-constrained situations. The second is a threshold to signal a switch from the Asp iterations to a point-source-only deconvolution for the rest of the deconvolution, meant to improve speed at a stage where Gaussian fits are no longer beneficial and a point-source model will suffice.
A test dataset was constructed with characteristics that are sometimes a challenge for CLEAN-based algorithms, but for which the ASP algorithm is expected to be more robust towards. The Hogbom, Multiscale and Asp algorithms were run, and metrics evaluated against a known true sky model. The control parameters of the Asp implementation were also exercised in order to ascertain whether they resulted in the expected changes in convergence behaviour.
Below is a summary. Please navigate using the links below to see details and plots. The "Comparison Results" links below contain the most concise information.
Observation Setup : VLA D-config at 1-2GHz with 5 channels and integrations spanning -5h to +55h Hourangle. Noise has been added.
Sky Model : A core-jet structure with spatial scales ranging from a point source to extended emission matched to the size of the uv-hole, and smooth continuum spectral structure ranging from a spectral index of 0.0 at the core, to about -1.0 in the extended lobe structure. A point source is also located atop the extended emission to simulate what is usually a difficult imaging case.
Numerical metrics include residual RMS, Chi-square and image-fidelity applied (as appropriate) to the residual image, restored image, and fitted spectral index image, and compared with a known truth image. Operational metrics include total runtime, the number of iterations done, and the number of major cycles triggered. All image cubes (including the truth model image used for the simulations) are smoothed to the same target resultion before comparison.
Test Hogbom/Multiscale/Asp with masks
Test Hogbom/Multiscale/Asp with no masks
Future Work/Tests : Additional topics/features to add/test.
Note : These results are based on one simulated dataset and so generalization of the results must be done with appropriate caution.
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
#%matplotlib widget
#%matplotlib inline
makeMSFrame(msn='sim_data', tel='VLA')
XPlot(vis='sim_data_VLA.ms', ptype='plotants')#,forceconvert=True)
XPlot(vis='sim_data_VLA.ms', ptype='uvcov')
listobs_jupyter(vis='sim_data_VLA.ms')
Two other spatial structures are available for these tests (in the supporting script). Note that the round mask size may need to be edited for these other examples.
Source structure - "basic" Three well-separated sources, two of which a point sources and one a 7arcmin Gaussian. This test the two extremes of angular resolutions.
Source structure - "simple" Two control point sources and a pair of blended Gaussians, each 5arcmin in size and with opposite spectral indices. This example is more relevent for a wideband reconstruction than a cube imaging exercise because it contains a turnover spectrum in the overlap region of the two Gaussians.
#source = 'basic'
#source = 'simple'
source = 'jet'
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
## Make the component list
makeCompList(clname_true='sim_'+source+'_vla.cl',tel='VLA',stype=source)
## Make an empty CASA image
makeEmptyImage(imname_true='sim_'+source+'_vla.im',tel='VLA')
## Evaluate the component list onto the CASA image
evalCompList(clname='sim_'+source+'_vla.cl', imname='sim_'+source+'_vla.im')
## Smooth to an angular resolution of 100 arcsec.
## This matches the Cube lowest resolution (first channel)
#smoothModel(imname='sim_'+source+'_vla.im',rbeam='100arcsec')
imsmooth(imagename='sim_'+source+'_vla.im',
outfile='sim_'+source+'.cube',
major='100.0arcsec',minor='100.0arcsec',pa='0.0deg',
targetres=True,overwrite=True)
fit_spectrum(cubename='sim_'+source+'.cube',
intensity='sim_'+source+'.cont',
alpha='sim_'+source+'.alpha',
beta='sim_'+source+'.beta',
pixt=0.1)
These average intensity, spectral index and curvature maps are to be used as a reference, to evaluate the quality of Cube (spectral line) imaging of the ASP (and MSClean and Hogbom) algorithms. They represent the best-case scenario of being able to reconstruct spectral index and curvature just from Cube imaging (and without MTMFS).
They also represent a low-resolution version of what ought to be deliverable from the current MTMFS implementation and the future Wideband ASP algorithm.
Simulate visibilities for the true sky model, using the standard gridder (no primary beams)
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
predictImager(msname='sim_data_VLA.ms',
imname_true='sim_'+source+'_vla.im',
tel='VLA',
gridder='standard')
copyModelToData(msname='sim_data_VLA.ms')
XPlot(vis='sim_data_VLA.ms', ptype='amp-freq-uvdist') #,forceconvert=True)
### Add noise.... later
addNoiseSim(msname='sim_data_VLA.ms',noise=1.0)
The following metrics were evaluated and then visualized to test several expectations about the behaviour of the Asp algorithm.
Let the a smoothed true image be $M$, output restored image be $I$, and output residual image be $R$
The RMS of the residual image pixel values : $rms(R)$
A goodness-of-fit metric evaluated against the smoothed true image : $\sqrt{\frac{\sum_{i=1}^{N} |I_i - M_i|^2}{N}}$ .
This metric is sensitive to differences across the entire image, giving equal weight to all pixels. A perfect fit will yield a value that matches the image rms noise level. We set $W=1$ over the entire flux density image.
An image fidelity metric defined as $1 - \frac{ \sum_{i=1}^{N} B_i W_i |M_i - I_i|}{\sum_{i=1}^N B_i W_i}$
where $W$ is a window function to select pixels of interest (e.g. PB cutoff mask), and $B$ is a pixel weight calculated as $B_i = max(|I_i|, |M_i|)$. The weight image $B$ emphasizes differences from the truth image for on-source pixels as well as away from the source. Such a weighted average of the absolute difference between the truth and output images is sensitive to off-source imaging artifacts. We set $W=1$ over the entire flux density image. (This is a modified version of eqn 10 from ngVLA memo 67).
The above metrics of chi-square and image fidelity are computed for the fitted spectral index from the output cube, using a window $W$ that encompasses only on-source pixels for which the spectral fit was itself successful. With complex extended emission, a dominant source of error in spectral index fits on cubes is the deconvolution inaccuracy per channel. These metrics assess the consistency of the reconstructions across spectral channels and its effect on a spectral index fit.
We record convergence history of the peak residual and total model flux per channel, as well as the total number of iterations done, the number of major cycles triggered, and the total runtime. The metrics provide a way to analyse the speed of convergence (numerical accuracy of each step) and associated compute cost.
For each cube imaging run, we display the reconstructed and residual images, convergence plots (per channel), the spectral index that may be fit from the output cube, and the relative difference between the smoothed true sky model and the reconstruction (normalized by the peak of the smoothed true image).
## A flag to run the notebook without re-running the imaging, but to only recompute metrics
run_imaging=False
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
all_results={}
if run_imaging==True:
imname='try_cube_hogbom'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='hogbom',niter=100000, gain=0.2, threshold='0.05Jy',nsigma=1.5,
interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
### Hogbom will produce better results if allowed to clean 'deeper' but then it will clean into
### the noise all around the source too. It needs automasking to reach on-source thresholds that
### ms-clean and asp achieve with out the auto-mask.
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
all_results['try_cube_hogbom'] = inspectResults(imagename='try_cube_hogbom',truthname='sim_'+source,maskthreshold=0.5)
if run_imaging==True:
imname='try_cube_multiscale'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='multiscale',niter=20000,gain=0.2, threshold='0.05Jy',nsigma=1.5,
scales=[0,3,10,20],interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
all_results['try_cube_multiscale'] =inspectResults(imagename='try_cube_multiscale',truthname='sim_'+source,maskthreshold=0.5)
if run_imaging==True:
imname='try_cube_asp'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='asp',niter=20000,gain=0.8, threshold='0.05Jy',nsigma=1.5,fusedthreshold=0.05,
largestscale=10, interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
all_results['try_cube_asp'] =inspectResults(imagename='try_cube_asp',truthname='sim_'+source,maskthreshold=0.5)
This test is to provide a comparison with the current usage mode of the ALMA pipeline.
if run_imaging==True:
imname='try_cube_hogbom_automask'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='hogbom',niter=100000, gain=0.2, threshold='0.05Jy',nsigma=1.5,
interactive=0,usemask='auto-multithresh',noisethreshold=4.0,sidelobethreshold=1.0,pbmask=0.01)
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
all_results['try_cube_hogbom_automask'] = inspectResults(imagename='try_cube_hogbom_automask',truthname='sim_'+source,maskthreshold=0.5)
print(all_results.keys())
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
compareResults(all_results)
These results show a clear improvement in numerical accuracy from Hogbom to MultiScale to Asp. This is true for the flux density as well as the fitted spectral index (alpha). The faster convergence of the Asp algorithm is also seen from the reduction in the number of minor cycle iterations and major cycles needed to reach the desired nsigma convergence criterion. This shows that the numerical behaviour for a working run, is exactly as expected.
The runtime of Asp is significantly higher than Hogbom and Multiscale. This too is expected since each iteration itself involves a fitting process to obtain best-fit Gaussian parameters. The relative runtime difference with Hogbom and Multiscale is only an illustration, derived from this one dataset.
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
nomask_results={}
if run_imaging==True:
imname = 'try_cube_hogbom_nomask'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
## 50k iterations to do the best it can
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='hogbom',niter=100000, gain=0.2, threshold='0.05Jy',nsigma=1.5,
interactive=0)
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
nomask_results['try_cube_hogbom_nomask'] = inspectResults(imagename='try_cube_hogbom_nomask',truthname='sim_'+source,maskthreshold=0.5)
if run_imaging==True:
imname = 'try_cube_multiscale_nomask'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='multiscale',niter=20000,gain=0.2, threshold='0.05Jy',nsigma=1.5,
scales=[0,3,10],interactive=0) ## Need to leave out the larger scale of 20 to avoid divergence.
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
nomask_results['try_cube_multiscale_nomask'] =inspectResults(imagename='try_cube_multiscale_nomask',truthname='sim_'+source,maskthreshold=0.8)
if run_imaging==True:
imname = 'try_cube_asp_nomask'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='asp',niter=20000,gain=0.2, threshold='0.05Jy',
nsigma=1.5,fusedthreshold=0.05,
largestscale=10, interactive=0)
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
## Asp displayed less tolerance of a high loop gain, without masks.
## gain=0.8 and cyclefactor=0.7 produced errors in the model image.
## gain=0.4 and cyclefactor was more stable.
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
nomask_results['try_cube_asp_nomask'] =inspectResults(imagename='try_cube_asp_nomask',truthname='sim_'+source,maskthreshold=0.5)
print(all_results.keys())
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
compareResults(nomask_results)
These results show that all three algorithms have a noticeable reduction in accuracy, without the stabilizing round mask. However, for one channel, the Asp algorithm had instabilities it could not recover from.
Initially, the Asp algorithm was tried with a loop gain of 0.8 but that run produced several spurious mid-to-large scale components in the model image in areas away from the source. Reducing the loop gain from 0.8 to 0.4 and then 0.2 (to match the setting for the other two algorithms) helped considerably but still did not generate adequate results for at least one channel. This suggests that (for this one example that is prone to instability) a simple region-based mask around the source was necessary for Asp.
Note : So far, this behaviour of Asp has been specific to this dataset.
Below is a comparison of the default behaviour of the implemented algorithm as well as the effect of using its two tuning parameters.
largestscale : By default, the implementation picks an initial guess of scale sizes that later evolve. largestscale limits the largest scale allowed for this initial guess. It should be picked as the largest scale for which appropriate data constraints may exist. It may be set, based of the uv-coverage of the observation (and not the actual sky structure as is required by the scale parameter in the msclean algorithm).
fusedthreshold : A threshold (on residual image pixel amplitude) below which the algorithm will use only point-source flux components. This is to reduce compute cost as well as to use delta functions when the residual pattern is no longer easily modeled by Gaussians.
tune_results={}
if run_imaging==True:
imname = 'try_cube_asp_default'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='asp',niter=20000,gain=0.2,threshold='0.05Jy',nsigma=1.5,
interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
tune_results['try_cube_asp_default'] =inspectResults(imagename='try_cube_asp_default',truthname='sim_'+source,maskthreshold=1.5)
## With largest scale set
if run_imaging==True:
imname = 'try_cube_asp_ls'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='asp',niter=20000,gain=0.2,threshold='0.05Jy',nsigma=1.5,largestscale=10,
interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
tune_results['try_cube_asp_ls'] =inspectResults(imagename='try_cube_asp_ls',truthname='sim_'+source,maskthreshold=0.5)
### With largest scale and a high loop gain
if run_imaging==True:
### Now try with largestscale set, and a high loop gain.....
imname = 'try_cube_asp_gain'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='asp',niter=20000,gain=0.8,threshold='0.05Jy',nsigma=1.5,largestscale=10,
interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
tune_results['try_cube_asp_gain'] =inspectResults(imagename='try_cube_asp_gain',truthname='sim_'+source,maskthreshold=0.5)
### With largest scale and a high loop gain and fusedthreshold set to something higher than the noise level...
if run_imaging==True:
imname = 'try_cube_asp_ft'
os.system('rm -rf '+imname+'.*')
casalog.setlogfile(imname+'.log')
st = time.time()
rec = tclean(vis='sim_data_VLA.ms',imagename=imname,imsize=512,cell='12.0arcsec',
specmode='cube',interpolation='nearest',nchan=5,start='1.0GHz',width='0.2GHz',
pblimit=-1e-05,deconvolver='asp',niter=20000,gain=0.8,threshold='0.05Jy',nsigma=1.5,largestscale=10,
fusedthreshold = 0.2, interactive=0,mask='circle[[256pix,290pix],140pix]')
et = time.time()
rec['runtime']=et-st
np.save(imname+'.summary.npy', rec)
tune_results['try_cube_asp_ft'] =inspectResults(imagename='try_cube_asp_ft',truthname='sim_'+source,maskthreshold=0.5)
%run -i Sim_Multiscale_Wideband_Script.py
%run -i Display_Experiments_Script.py
compareResults(tune_results)
Blue : With the defaults, for this dataset, Asp performs a very large number of iterations without much benefit (i.e. it did not reach minor cycle stopping thresholds as soon as the later tuned runs). The initial guess of scale sizes includes two sizes larger than 10pix, and could have contributed to unconstrained iterations. (Note : For this dataset, the Multiscale algorithm also exhibits a similar instability when provided with scale sizes beyond 10pix. Tests on other simulated and real datasets (not reported here) did not show such instabilities.)
Orange : With a limit on the largest scale, the convergence stabilized very well, converging faster and with better accuracy. This demonstrates the efficacy of this control parameter and that it is behaving as intended.
Green : The addition of a high loop gain (from the default of 0.1 to 0.8) showed faster convergence and also numerical improvement. This demonstrates the expected behaviour for an algorithm with more accuracy per minor cycle iteration.
Red : Switching to point source deconvolution partway through the deconvolution improved runtime (by 20% for this example and setup) but degraded the numerical performance only very slightly (in comparison to the other Asp runs). This demonstrates the expected behaviour of improved speed at the cost of some accuracy, of switching to a Hogbom-like algorithm below a certain threshold.